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Abstract 

Using Brownian dynamics simulations, we study the migration of long charged 
chains in an electrophoretic microchannel device consisting of an array of micro- 
scopic entropic traps with alternating deep regions and narrow constrictions. Such 
a device has been designed and fabricated recently by Han et al. for the separation 
of DNA molecules (Science, 2000). Our simulation reproduces the experimental ob- 
servation that the mobility increases with the length of the DNA. A detailed data 
analysis allows to identify the reasons for this behavior. Two distinct mechanisms 
contribute to slowing down shorter chains. One has been described earlier by Han 
et al. : the chains are delayed at the entrance of the constriction and escape with a 
rate that increases with chain length. The other, actually dominating mechanism is 
here reported for the first time: Some chains diffuse out of their main path into the 
corners of the box, where they remain trapped for a long time. The probability that 
this happens increases with the diffusion constant, i. e., the inverse chain length. 

Key words: gel electrophoresis, microfluidic system, DNA separation, entropic 
trap, computer simulation 



1 Introduction 



Electrophoresis is the standard method of DNA separation by length (Viovy, 
2000). Since the mobility of DNA molecules in free solution is independent 
of their length, the DNA is commonly introduced into a gel, which serves as 
a random sieve. Unfortunately, the efficiency of gel electrophoresis decreases 
with the length of the DNA. Moreover, bioanalytic devices are more and more 
miniaturized, and incorporating random gels with characteristic pore sizes in 
the nanometer range into future nanodevices will presumably become increas- 
ingly problematic. Therefore, much recent effort has been devoted to designing 
and making well-defined microstructured devices for DNA separation (e. g., 
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Fig. 1. Schematic cartoon of the microchannel device proposed by Han et al. (2000) 
(side view). The width of the channel in the y-direction is much larger (~ 30/xm). 



Han et al, 1999; 2000; 2002; Bader et al., 1999; Hammond et al, 2000; Duong 
et al., 2003). 

In this paper, we focus on a geometry proposed by Han and coworkers (Han 
et al., 1999), which is based on the idea of entropic trapping. The DNA is in- 
troduced into a small microchannel with alternating deep regions and shallow 
constrictions. The channel thickness in the constrictions is much smaller than 
the radius of gyration of the DNA molecules. The deep region is large enough 
to accommodate the equilibrium shape of the DNA molecules. Entering the 
constrictions is thus entropically penalized, and the deep regions can be in- 
terpreted as entropic traps. A schematic cartoon of the structure is shown in 
Fig. 1. 

In these structures, Han et al. made the counterintuitive observation that long 
DNA molecules travel faster than short molecules. They explained their find- 
ing by means of a simple kinetic model. In order to escape from one of the deep 
regions, the DNA must overcome an activation barrier AF C . The escape pro- 
cess is initiated by a thermally activated stretching of a chain portion (length 
x) into the shallow constriction. This costs entropic penalty proportional to x, 
but it also decreases the electric potential energy by an amount ~ Ex 2 , where 
E is the electric field. Hence there exists a critical length x c oc 1/E, below 
which the entropy penalty dominates and the chain is driven back into the 
deep region. At x > x c , the energy gain dominates and the chain escapes. The 
free energy at x c represents the activation barrier for the escape process. It 
depends solely on the electric field, AF C oc 1/E. The rate of escape attempts, 
1/tq, increases with the chain size, since the amount of polymer in contact 
with the constriction is larger for larger chains. The mean trapping time in 
this simple model is given by 

r = r exp(AF c /k B T), (1) 



where T is the temperature and k B the Boltzmann factor. This implies that 
long chains travel faster. 

The results of Han et al. have motivated recent Monte Carlo simulations by 
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Tessier et al. (2002). These authors used the Bond fluctuation model, a well- 
known lattice model for polymers, to study the mobility of polymers in a 
geometry similar to that suggested by Han et al. within local Monte Carlo 
dynamics (single monomer moves). Their results confirm the trapping picture 
of Han et al., and even details of the kinetic model. For example, they present 
evidence that the penetration depth x of the chain into the constriction can be 
used as a "reaction coordinate" for escaping, with a critical value x c oc 1/E. 

When looking more closely at the simulation data, one notices that the trap- 
ping in the system of Tessier et al. is unexpectedly strong. Although the width 
of the constriction is more than twice as large as in the experiments (in units 
of the persistence length), the molecules spend almost the entire time in a 
trapped state, and very little time traveling, even at intermediate fields. This 
effect might be an artifact of the lattice model. We note that the width of 
the constriction in the simulations is just 10 lattice constants, while every 
monomer occupies a cube with 8 lattice sites, and the average bond length 
is 2.8 lattice constants. Moreover, the Monte Carlo dynamics is not realistic. 
Monomers are picked and moved randomly (with some moves being rejected), 
whereas in the real system, they are pulled into the constriction by the electric 
force. Dynamic Monte Carlo is known to reproduce diffusional and relaxational 
dynamics in systems near local equilibrium. Nevertheless, it is not clear how 
well it can be applied to study driven systems. The details of the dynamics 
should matter particularly at higher electric fields, in situations where the 
chains travel so fast that they no longer reach local equilibrium in the traps. 
Therefore, simulations of off-lattice models with a more realistic dynamical 
evolution are clearly desirable. 

As a first step in this direction, Chen and Escobedo (2003) have recently 
studied the free energy landscape of chains in a single, non-periodic trap with 
Monte Carlo simulations of an (off-lattice) bead-spring model. The initial con- 
figuration in these simulations is that of a fully relaxed chain in the absence of 
an electric field. With that starting point, the free energy barrier AF turns out 
to depend on the chain length for short chains, and to level off at higher chain 
lengths. The data for AF do not seem to support the relation AF oc 1/E. 
Since the simulations still used Monte Carlo, the results on dynamical prop- 
erties were limited. 

In the present paper, we present Brownian dynamics simulations of a full 
(non-equilibrium) periodic array of entropic traps. We employ a Rouse-chain 
model, such as has been successfully used by others to study the migration of 
DNA in various geometries (Deutsch, 1987; 1988; Matsumoto, 1994; Noguchi, 
2001). Our main result is the observation of a new trapping mechanism in 
geometries such as Fig. 1, which is at least as important as that suggested by 
Han et al., and which can also be exploited to achieve molecular separation. 
The insight gained from our study should be useful for the design of future, 
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improved molecular separation devices. 



The paper is organized as follows: in the next section, we describe and discuss 
the model and give some technical details on the simulation. The results are 
presented in section 3. We summarize and conclude in section 4. 



2 The Simulation Model 

We model a single DNA molecule as a chain with pairwise interactions 

{a/r) 12 -(a/rf + l/A : (r/af < 2 



V P air(r)/k B T = < 



(2) 

: otherwise. 



This potential describes soft, purely repulsive interactions between beads of 
diameter a. The beads are connected by springs with the spring potential 

k 2 

V S pring(r) = —T . (3) 



The spring constant was chosen very large, k = lOO/c^T/a 2 , in order to prevent 
chain crossings. 

The molecule is confined in a structured channel with a geometry similar to 
that used by Han et al.. The thickness of the shallow constrictions and thick 
regions is t s = 7a and t t = 80a, respectively, the length of a deep region is 80 
a, and the total length of a period is L = 100a. In the lateral direction, the 
channel is infinite. The chain beads are repelled from the walls of the channel 
by means of a wall potential essentially identical to (2), Vy/aiiij') = V pa i r (r), 
where r is the distance of a bead to the closest wall point 1 . The details of the 
wall potential do not influence the results, as long as it is repulsive and short 
ranged. Note that the effective width of the channel, i. e., the width of the 
space accessible to a bead, is reduced by 2a with this potential. A snapshot 
of a chain in such a channel is shown in Fig. 2. 

DNA is a charged polyelectrolyte with 2e~ per base pair, thus each bead carries 
a charge q and is subject to an electric field E = — V$. We assume that the 
charges do not interact with one another. This will be discussed further below. 
The distribution of the electric potential $(r) in the channel was calculated 
numerically by solving the Laplace equation (A$ = 0) with von Neumann 
boundary conditions at the walls (n ■ V$ = 0, where n is the surface normal). 



At corners, the potentials of the adjacent walls are summed up 
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Fig. 2. Snapshot of a chain with N = 1000 beads in an entropic trap. The solid 
lines show the electric field lines. 



The dynamical evolution of the system is described by a Kramer's equations 
(Risken, 1989) 

*i = v, (4) 
mv, = fj - Cvj + rji, 

where fj is the total force acting on bead i, and Vj its velocity, and ( is a 
friction coefficient. The random noise r\ fulfills 

(rha) = Q (5) 
(Via(t)Vjp(t')) = 2k B TCS ij 5 al3 5(t - f), 

with a, (3 — x, y or z. The random noise and the friction mimics the effect 
of the solvent surrounding the DNA. The chains differ from standard Rouse 
chains (Doi, 1986) by their excluded volume interactions and by the inertia 
m of the beads. With the parameters of our model (m = l.(, 2 a 2 /k B T , see the 
discussion further below), the chains exhibit Rouse dynamics on length scales 
of the order the gyration radius of the chain, R g (Streek, 2002). 

Several important physical factors are neglected in the model. First, electro- 
static interactions within DNA chains are not taken into account. This is 
justified by the fact that the Debye screening length in typical electrophoresis 
buffers is only a few nm, comparable to the diameter of the double helix. Sec- 
ond, we do not consider electro-osmotic flow, or flow in general. This reflects 
the experimental situation reported by Han et al.. Third, hydrodynamic inter- 
actions are disregarded. This is partly justified by the theoretical observation 



that for polyelectrolytes dragged by an electric field, the hydrodynamic inter- 
actions should be screened over distances larger than the Debye length, since 
the counterions are dragged in the opposite direction (Viovy, 2000). Unfortu- 
nately, the argument is only strictly valid as long as no non-electric forces are 
present (Long, 1996). In our case, the walls of the channels exert non-electrical 
forces which stop the polymers, but do not prevent the counterions in the De- 
bye layer from moving. Thus the chains experience an additional trapping 
force from the friction of the counterions. This effect is disregarded in our 
model. Furthermore, diffusion is not treated correctly even in free solution. 
The diffusion constant of Rouse chains scales as 1/N with the chain length N. 
Including hydrodynamic interactions, one expects Zimm scaling 1/R g , where 
the gyration radius R g scales like R g oc iV 3 / 5 for self-avoiding chains. Experi- 
mentally, the diffusion constant of DNA in typical buffer solutions is found to 
scale as l/iV a672 (Stellwagen, 2003). 

Unfortunately, a full treatment which accounts correctly both for the (dynami- 
cally varying) counterion distribution as well as the hydrodynamic interactions 
is not feasible with standard computational resources. The simplifications of 
our model influence the results quantitatively, but do presumably not change 
them by orders of magnitude. The qualitative behavior should not be affected. 

Keeping these caveats in mind, we can now proceed to establish a quantitative 
connection with the experimental setup of Han et al.. The natural units in our 
model are related to the parameters a, (, \q\ (the charge per bead), and T (the 
temperature). More specifically, the energy unit is e = ksT, the length unit is 
a, the time unit is to = (a 2 /ksT, and the electric field unit is E = ksT ja\q\. 
Throughout the paper, all quantities shall be given in terms of these natural 
units. They shall now be related to real (SI) units. 

Since the experiments are carried out at room temperature, the energy unit 
is e = 300k bK = 4 • 10~ 21 J. The length unit is obtained from matching the 
persistence length of the model chains, l p = 1.6a (Streek, 2002) with that of 
DNA, l p = 45nm, yielding a = 30nm. The persistence of the chain is also used 
to determine the number of base pairs (bp) per bead. A DNA molecule contains 
approximately 150 bp on a stretch of length l p . In our model, the average bond 
length between two beads is 0.847a (Streek, 2002), thus we have 1.9 beads per 
persistence length, and one bead represents roughly 80 bp. The elongation for 
a chain crossing with minimal energy is r cross as 1.198(7, which corresponds 
to an energy barrier AE CIOSB ss 82E for crossing. The Boltzmann-factor for 
this energy barrier turns out to be less than 10~ 35 . To check the simulation 
program, the bond length distribution was compared to the Boltzmann-factor 
and very good agreement was found. Furthermore, no bond ever exceeded a 
length of l.lcr. Thus no chain crossing occured in our simulations. 
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The time scale to is calculated from the diffusion constant D. For Rouse chains 
of length N, D is given by 



D = 



k B T a 2 



(6) 



iVC Nt ' 



Experimentally, Stellwagen et al. (2003) have recently reported the relation 



Choosing as reference a chain of length 40 kbp (N = 500), we obtain t = 
2.9 • 10~ 6 s. 

Finally, the unit of the electric field can be identified from matching the mo- 
bility fiQ of free chains. The theoretical value is fio = \q\/C = °"Ao-^o- 111 
experiments, the mobility depends strongly on the choice of the buffer. Un- 
fortunately, Han et al. do not report explicit measurements of the free-chain 
mobility of DNA in the buffer used in their experiments. In the microchan- 
nel, they observe that the overall mobility saturates at high field strengths. 
The apparent maximum mobility /i max results from an average over a slow 
motion in the reduced electric field of the deep regions, and a fast motion in 
the enhanced electric field of the constrictions (see Fig. 2). For large periods 
L ^> ti,t s , the ratio of these two field strengths is simply the inverse of the 
thickness ratio a = U/t s . One can then derive a relation between the true free 
chain mobility /io and the apparent mobility fi max - 

A*o = /Wr (xi + x s /a) (xi + x s a) , (8) 

where xi and Xg £1X6 the relative lengths for the deep regions and shallow 
constrictions, respectively (xi + x s = 1). Han et al. measure fi ma x ~ 0.13 • 
10~ 8 m 2 /Vs in an experimental setup with x\ = x s = 1/2, a ~ 15, and L = 
10 — 40/im (Han, 1999). Using Eq. (8), one can thus estimate /xq = 4.3/i max = 
0.55 • 10~ 8 m 2 /Us. This value of /x seems very low compared to typical values 
in the literature. Stellwagen et al. have measured DNA mobilities in Tris- 
acetate buffers at various concentrations and found values between 2 — 4 • 
10 _8 m 2 /Us (Stellwagen, 2002). In 40 mM Tris-acetate buffer, they obtain 
Hq = 3.3 • 10 _4 cm 2 /Us (Stellwagen, 2003). Using the latter value as an order- 
of-magnitude estimate, we identify E ~ 3. • 10 3 U/cm. However, the results 
of Stellwagen et al. can probably not be applied here, because the mobility 
depends strongly on the buffer. The first estimate, /i = 0.55 • 10~ 8 m 2 /Us, 
leads to the identification E ~ 2. • 10 4 U/cm. 

Han et al. (1999, 2000) have separated DNA of lengths between 5 and 160 kbp. 
Here we study chains of length iV < 1000, corresponding to < 80 kbp, which 



D = 7.73 ■ 10 6 (number of base pairs) 



-0.672 



cm s 



(7) 
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is comparable. The depth of the deep channel regions in the experiments was 
ti — 1 — 3/im, which also compares well with the simulation, 80a = 1.6/im. 
The depth of the constriction was t s = lOOnm 2l p in the experiments. In 
our case, it is 5a ~ 3l p , which is slightly larger, but still comparable. Since 
our channels are wider, we expect that the trapping in the simulations will 
be less pronounced than in reality. The total length per trap (period) was 
L = 4 — 40/im in the experiments (mostly 4/xm) and L = lOOo" = 3/im in the 
simulations. 

The average electric field strength in the experiments was varied between 20 
and 80 V/cm, which corresponds to ~ 0.001 - 0.004£ or ~ 0.006 - 0.03£ , 
depending on the identification of Eq. In the simulations, we studied field 
strengths between 0.0025 — 0.04i?o- When comparing field strengths, we must 
keep in mind that the electric field in the channel is not homogeneous (see 
Fig. 2). In our geometry, the electric field in the constriction is enhanced by a 
factor of 2.5. In the experimental setup, the enhancement factor is only ~ 1.8, 
due to the fact that the length ratio between the shallow and deep regions 
is different (50:50 in the experiment, 20:80 in the simulations). The ratio of 
field strengths in the shallow and deep regions in our simulations is roughly 
4. In the experimental geometry, it is ti/t s 10 — 30 for large periods L and 
smaller otherwise. 

The remaining model parameter that has yet to be determined is the mass m 
of a bead. We note that the actual value of the mass has no influence on the 
static properties of the chain (e. g., the chain flexibility), nor on the diffusive 
part of the dynamics. It does, however, determine the relative importance of 
vibrational modes in the chain and other inertia effects. The latter can be 
characterized by the electrophoretic relaxation time r e , i. e., the characteristic 
decay time of the drift velocity of free flow DNA, if the electric field is suddenly 
turned off. Typical values of r e are r e ps 10~ 9 — 10~ 12 s (Grossmann, 1992). In 
our model, the electrophoretic relaxation time is r e = m/(, and the correct 
value of the mass would be m ~ 10~ 3 — 10~ 6 £to- However, this is unfortunate 
from a computational point of view, because the mean velocity per bead, 



V = v ^ksT/m, becomes large for such small bead masses, and the time 
step has to be chosen short as a consequence. The simulation becomes very 
inefficient. On the other hand, we are not interested in inertia effects here, 
and we wish to study dynamical effects on much longer time scales than to- 
Therefore, we chose to give the beads an unphysically high mass, m = l^to, 
leading to an electrophoretic relaxation time r e = t ~ 10~ 6 s. On time scales 
t and less, the dynamics will thus be unrealistic, but this does not affect the 
slow diffusive dynamics. 

We close this section with a few technical remarks. The dynamic equations (4) 
were integrated using the Verlet algorithm. The stochastic noise was realized 
by picking in every time step a vector 77 at random. Diinweg and Paul (1991) 
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have shown that the distribution of random numbers in such a procedure does 
not necessarily have to be Gaussian. Here, we used a uniform distribution in 
a sphere. Since we consider single chains only, no periodic boundaries were 
necessary. With the mass m = l(t , the time step could be chosen 0.01 to- 
Typical run lengths were between 4 and 20 million t (5-25 seconds). The sim- 
ulation jobs were managed by the Condor Software Program (Condor), which 
was developed by the Condor Team at the Computer Science Department of 
the University of Wisconsin (Condor, 2003). 



3 Results and Discussion 

Three examples of trajectories at the average field E = 0.005Eq are shown 
in Fig. 3. They reveal three qualitatively different modes of migration. The 
trajectory of very short chains (N = 10 corresponding to 800 bp) is dominated 
by diffusion. The chains wander back and forth in the trap, until they even- 
tually escape into the next trap. The movement of chains with intermediate 
chain length (N = 100 or 8 kbp) is much more directional, but still irregu- 
lar. They are often trapped at the entrance of constrictions. In contrast, long 
chains (N = 1000 or 80 kbp) travel smoothly. The first (leading) monomer 
is sometimes trapped, but this does not arrest the rest of the chain. Whereas 
the time spent in one box fluctuates strongly for shorter chains, it is roughly 
constant for large chains. 

The behavior of intermediate chains (N = 100) has similarity to that observed 
in the simulations of Tessier et al. (2002). However, the trapping is much less 
pronounced in our case, and chains frequently pass from one box to another 
without being trapped at all. Trapping effects comparable to those reported 
by Tessier et al. were only observed at the lowest field strength, E = 0.0025-E - 
At that field value, chains of all lengths got trapped. As we shall see below, 
however, chain separation turned out to be not efficient for such small fields. 

We have carried out simulations for seven different chain lengths N = 10, 20, 
50, 100, 200, 500, and 1000, and for five average field values, E/E = 0.0025, 
0.005, 0.01, 0.02, and 0.04. The resulting mobilities, determined as /i — (v)/E, 
are summarized in Fig. 4. For all fields except the lowest, the mobility increases 
steadily with the chain length. At high chain lengths and high fields, it begins 
to saturate. The maximum value fi max is only about half as large as the free 
chain mobility /io, due to the fact that the chains spend a disproportionate 
amount of time in the deep regions, where the local electric field is lower than 
average. 

At the lowest field (E = 0.0025i?o), the mobility depends only slightly on 
the chain length and even decreases with chain length for small N. At such 



9 




Fig. 3. Trajectories of chains moving in an entropic trap array for three different 
chain lengths N in the presence of an average electric field E = 0.005£'o. The middle 
solid line shows the position of the center of mass. The (dashed) upper and lower 
lines show the positions of the first and the last monomer. In the case N = 10, 
these three lines cannot be distinguished from each other. The dashed horizontal 
lines indicate the positions where constrictions begin. 



small fields, backwards diffusion becomes important for short chains. We have 
seen in Fig. 3 that short chains explore the whole trap. At E — 0.0025-E > 
they sometimes even travel backwards into the trap that they just left. The 
conformational entropic penalty for entering the constrictions is small for short 
chains. In the limit E — > 0, the inverse mobility is thus simply proportional 
to the number of times the chain visits the entrance of the constriction. Since 
the latter scales like N (the inverse diffusion constant D^ 1 ), the mobility is 
then expected to decrease with chain length. In our system, this is observed 
at E = 0.0025^0 for chain lengths smaller than N = 100. For N > 100, the 
mobility increases again with chain length. The resulting overall chain length 
dependence is small, and the chain separation is not efficient. 

The quality of molecular separation systems is often characterized in terms of 
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Fig. 4. Mobility \x in units of the mobility /j,q of a free chain as a function of chain 
length N for different average electric fields E in units of Eq . The dashed lines show 
the prediction of Eq. (12) for the fields E/E = 0.04,0.02,0.01,0.005 (from top to 
bottom). 

the theoretical plate number 

Opiate = 16 (W%) 2 - (9) 



where tn is the retention time, i. e., the total time spent in the system, and 
t w the width of the peak at the baseline. Fig. 5 shows the plate number per 
trap for our system. In the interesting regime, we have 10-100 plates per trap. 
At ~ 10 5 traps per meter, this corresponds to theoretical plate numbers of 
10 6 — 10 7 plates/m, which is quite good and in agreement with the results of 
Han et al. (2000). 

We will now investigate the migration modes in more detail. To this end, we 
have calculated histograms of the retention times spent in one trap. They were 
defined as the difference t n+ \ — t n of the times t n when the first monomer of 
a chain first enters the deep region of the nth trap. Fig. 6 shows distributions 
of retention times for chains of different length N in the field E = 0.01-Eb- 

The chains need a minimum time ~ 1 — 1.5 • 10 4 to to travel from one trap to 
another. After that "dead" period, the histogram rises rapidly and reaches a 
maximum at t max ~ 2. • 10 4 to- iFrom a comparison of histograms for all our 
simulation data (not shown), we deduce that the position of the maximum 
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Fig. 5. Theoretical plate number per trap as a function of chain length for different 
electric fields E as indicated (in units Eg). 



tmax of the distribution depends very little on the chain length and is strictly 
proportional to the inverse field 1/E. The product L/t max E = fi max ~ 0.5/zo 
determines the maximum mobility fi m ax m our system. 

Beyond the maximum, the distribution decays rapidly for long chains (N = 
1000), and much more slowly for short chains (N = 10), following an exponen- 
tial behavior. This is consistent with the common picture, where the migration 
is determined by one single escape rate 1/r. The histogram for intermediate 
chain length iV = 100, however, reveals that the situation is more complex in 
reality. The initial decay of the distribution can be fitted with one exponen- 
tial, however, a long-time tail emerges at times beyond t ~ 5 ■ 10 4 to- Thus the 
distribution of retention times at N = 100 has two characteristic time scales 
Tfast and t s1ow . 

To some extent, this phenomenon is already apparent in the trajectory of 
Fig. 3. The figure suggests that there exist two qualitatively different ways 
how chains pass from one trap to another: Either they travel relatively straight 
and unimpeded, or they get trapped and linger for some time at the border 
of the constriction. 

The trapping mechanism suggested by Han et al. alone cannot explain these 
observations. Here we propose an additional trapping mechanism, which also 
slows down short chains and which presumably accounts to a large extent for 
the chain length dependence of the mobility observed in our simulations. The 
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Fig. 6. Distribution of Retention times in one trap for different chain lengths N 
and electric field E = O.OlE'o. The thick solid line is a fit to the initial exponential 
decay at chain length N = 100. 



idea is that chains get trapped at the side walls and corners of the deep boxes 
due to diffusion. The electric field lines lead the chains from the outlet of one 
constriction directly to the entrance of the next one. With a certain probability, 
the chain will therefore reach the next constriction without detours, and then 
get delayed there due to the mechanism suggested by Han et al.. This accounts 
for the fast time scale r fast . On the other hand, chains may also diffuse out of 
the main path. They may access regions of the trap at the walls and in the 
corners where the electric field is very small. In that case, they are caught in 
a force-free region, and they can "escape" only by diffusion. 

We will now explore whether our data support this picture. In the following 
analysis, only data for chain lengths iV > 20 and fields E > 0.005E were used. 
In most of these systems two time scales r were observed. The value of the fast 
time scale r fast could be extracted by fitting an exponential (Aexp(— t/r fast )) 
to the initial decay of the histogram H(t). The determination of the slow 
time scale was much more difficult, due to the poor statistics for the late 
time tails of the histograms. We used two approaches: First we fitted the long 
range tail with an exponential function to obtain a rough estimate. Assuming 
H(t) oc e"'/ rslow , we calculated r slow via 

oo oo 

r slow = J dt(t-t cut )H(t) I J dtH(t) (10) 
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Fig. 7. Characteristic length scales (rescaled) of the retention time distribution 
as a function of chain lengths N for different electric fields E. The filled symbols 
correspond to the fast time scale r fast , the open symbols to the slow time scale r slow . 
The thick lines show for comparison power laws as indicated. 



which is independent of t cut . We used t cut = 500/ E (to-E'o) to analyze the data. 
We checked whether the result deviated strongly from the previous estimate, 
and whether it depended strongly on the cutoff t cut . If this was the case, 
(because the data for H(i) scattered strongly), the result was discarded. The 
remaining values are compiled in Fig. 7, together with the data for r fast . 

If our suspicion is correct that the fast process corresponds to the mechanism 
of Han et al., then the rate l/r fast should be proportional to the amount of 
polymer in contact with the channel. Since the channel entrance is essentially 
one dimensional, the contact area should be proportional to the gyration radius 
R g of the chain. Fig. 7 shows that the fast time scale indeed scales like r fast oc 
l/R g oc N- 3 / 5 . 

The relation between r fast and the electric field strength E is more complicated. 
According to Han et al., the chains must overcome a free energy barrier of 
height proportional to 1/E in order to escape. On the other hand, the prefactor 
To in Eq. (1) may also depend on E. The resulting ^-dependence can be rather 
complex. In the field range of our simulation, the resulting relation r fast (E) can 
be approximated by the empirical law r fast oc E~ 1,55 . 

Like the fast time scale r fast , the slow time scale r slow also decreases with the 
chain length, but the dependence here is very weak. Unfortunately, the quality 
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Fig. 8. Schematic cartoon of a chain caught at the corner of the trap in the field- free 
region. The chain experiences a net force towards the wall. 

of the data was not sufficient for a quantitative analysis. If our interpretation 
is correct, then r slow characterizes an escape from a low- force region. We note 
that extended chains in a corner experience a net electric force towards the 
wall, even though the field lines of course never enter the wall (see Fig. 8). 
In order to leave the corner, the chain must either move against the force, 
or change its shape. In both cases, it has to overcome a free energy barrier 
before it can get rescued from the field lines. A simple Ansatz would predict 
the escape probability to be proportional to the diffusion constant D oc l/N, 
and to the area covered by the chain, R 2 oc iV 6 / 5 , or the total chain length N. 
This would yield a very weak net chain length dependence r slow ~ iV^ 1 / 5 or 
Tsiow ~ N°, which is consistent with the observed behavior (Fig. 7). 

Thus the slow time scale r slow itself does not contribute much to the chain 
length dependence of the mobility. The main effect comes from the fact that 
the relative number of chains caught in the field-free region depends on the 
chain length. The travel time from one channel to the next for undeflected 
chains is slightly smaller than t max = L/E/i max . We assume that chains have 
to diffuse at least over a distance z into the field-free region (direction —z, 
see Fig. 1) in order to get caught. The distribution along the ^-direction after 
a time t max will be approximately Gaussian: N(z) oc e~ z ^ 6Dt . The probability 
of being caught is thus 

oo 

P = J N(z) = etfc{z /^6Dt max ) = erfc(aVNE), (11) 

20 

where erfc(y) = dxexp(— x 2 ) is the complementary error function. 

We can test this prediction under the assumption that the two time scales r fast 
and r slow are sufficiently far apart that they can be separated. In that case, 
one can choose a time t cut such that almost all undeflected chains have left the 
trap, while almost all deflected chains are still in the field-free region. P can 
then be approximated by the relative number P of chains left in the field-free 
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Fig. 9. Probability Po(t) that a chain is still in the trap at the time 
i cut = 350/E (toEo), vs. NE in units of Eq. The solid line represents a fit to Eq. (11) 
with a = 0.59. 



region at the time i cut . Fig. 9 shows the result of such an estimate. The data 
collapse reasonably well for different N,E and can be fitted with Eq. (11), 
with the fit parameter a = 0.59. Inserting fi max = 0.5/i and L = 100a, one 
obtains z ~ 20a. Thus chains get caught if they sidetrack by more than ~ 20a 
from their main path, which is determined by the field lines. 

These considerations establish the existence of a new mechanism which pro- 
duces chain length dependent mobility. To assess the relative importance of 
the new mechanism, we compare the real mobility data, Fig. (4), with a very 
simple model. Chains either travel straight across the trap, or get sidetracked 
into the field-free region. Traveling across the trap takes at least the time 
tmax = LE/fi max . The chains caught in the field- free region spend an ad- 
ditional time r slow in the trap. We make the simplification that Et s1o „ is in- 
dependent of the chain length and electric field and given by the number 
Er slo „ = 400 t E Q , which is roughly the value at chain length ~ 100 and field 
strength ~ 0.01 — 0.02 Eq. The relative number of chains caught in the deep 
region P is calculated according to Eq. (11), with a = 0.59 (taken from Fig. 9). 
The resulting mobility is 

(12) 



A*0 Mr 



+ etfc(aVNE) 
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This prediction is compared with the actual data in Fig. 4 (dashed lines). 
Despite the simplicity of Eq. (12), the agreement is remarkably good. 



4 Summary and Outlook 

To summarize, we have presented the first off-lattice Brownian dynamics sim- 
ulation of DNA migration in an entropic trap array. We reproduce the ex- 
perimentally observed phenomenon that the mobility increases with the chain 
length. This result can be traced back to two distinct mechanisms. The first 
mechanism has already been discussed by Han et al.: Chains get delayed at 
the narrow channels connecting the traps. They escape through the channels 
with a probability which is proportional to the radius of gyration of the chain 
and thus scales as iV 3//5 . However, we found that this effect accounts only 
in part for the total chain length dependence. In the second mechanism, the 
chains are trapped with a certain probability at the side and in the corner 
of the box. The characteristic time for escaping such a configuration is very 
long. The trapping probability increases with the diffusion constant, which is 
in turn inversely proportional to the chain length. As a result, the mobility 
increases with the chain length. 

To our best knowledge, this mechanism has not yet been described in the 
literature. It becomes relevant when the period L of the structure is small. 
Indeed, Han et al. have studied structures with periods ranging from 4 to 40 
/im (Han, 1999), but they reported separation by length only for their smallest 
structure with L = 4/im, in a system with dimensions comparable to those 
studied here. 

We have observed a number of other phenomena, which we shall not describe 
in detail here. In contrast to Chen et al. (Chen, 2003), we have considered truly 
non-equilibrium systems. Subsequent escapes cannot necessarily be considered 
as independent events. The longest chains N = 1000 do not recover the equi- 
librium coil structure in the middle of the trap, but they remain stretched in 
the x direction. Moreover, our data seem to provide evidence that chains even 
retain some memory of the previous escape process. At high fields, successive 
escape times seem to be correlated. Unfortunately, the statistical quality of 
the data is not good enough to allow for a more thorough analysis. 

The situation becomes even more complicated when the shallow channel is 
made wider. Whereas in the cases presented here, long chains migrated faster 
than short chains, we have observed the inverse effect in microfluidic system 
with wider channels (Duong, 2003). The effect as well as other, even more 
unexpected phenomena, can be reproduced in simulations. These phenomena 
will be presented in a forthcoming publication (Streek, 2004). 



17 



Acknowledgments 

We thank Ralf Eichhorn for critically reading the manuscript. This work was 
funded by the German Science Foundation (SFB 613, Teilprojekt D2). 



References 

Bader, J. S., Hammond, R. W., Henck, S. A., Deem, M. W., McDermott, G. A., 
Bustillo, J. M., Simpson, J. W., Mulhern, G. T., Rothberg, J. M., 1999. 
DNA transport by a micromachined Brownian ratchet device. PNAS 96, 
13165-13169. 

Chen, Z., Escobedo, F. A., 2003. Simulation of chain-length partitioning in a 
microfabricated channel via entropic trapping. Mol. Sim. 29, 417-425. 

The Condor team, 2003. software package from www.cs.wisc.edu/condor/. 

Deutsch, J. M., 1987. Dynamics of pulsed-field electrophoresis. Phys. Rev. 
Lett. 59, 1255-1258, 

Deutsch, J. M., 1988. Theoretical studies of DNA during gel electrophoresis. 
Science 240, 922-924. 

Doi, M., Edwards, S. H., 1986. The theory of polymer dynamics. Clarendon 
press, Oxford. 

Diinweg, B., Paul, W., 1991. Brownian dynamics simulation without Gaussian 
random numbers. Int. J. Mod. Phys. C 2, 817-827. 

Duong, T. T., Kim G., Ros, R., Streek, M., Schmid, F., Brugger, J., Ansel- 
metti, D., Ros, A., 2003. Size-dependent free solution DNA electrophore- 
sis in structured microfluidic systems. Microelectronic Engineering 67-68, 
905-912. 

Grossmann, P. D., Colburn, J. C, 1992. Capillary Electrophoresis, Theory 
and Practice: Free Solution Capillary Electrophoresis. Acad. Press, San 
Diego. 

Hammond, R. W., Bader, J. S., Henck, S. A., Deem, M. W., McDermott, 
G. A., Bustillo, J. M., Rothberg, J. M., 2000. Differential transport of 
DNA by a rectified Brownian motion device. Electrophoresis 21, 74-80. 

Han, J., Turner, S. W., Craighead, H. G.. 1999. Entropic trapping and escape 
of long DNA molecules at submicron size constriction. Phys. Rev. Lett. 
83, 1688-1691; Erratum. 2001. Phys. Rev. Lett. 86, 1394. 

Han., J., Craighead, H. G., 2000. Separation of long DNA molecules in a 
microfabricated entropic trap array. Science 288, 1026-1029. 



18 



Han., J., Craighead, H. G., 2002. Characterization and optimization of an 
entropic trap for DNA separation. Anal. Chem. 74, 394-401. 

Long, D., Viovy, J.-L., Ajdari, A., 1996. Simultaneous action of electric fields 
and nonelectric forces on a polyelectrolye: motion and deformation. Phys. 
Rev. Lett. 76, 3858-3861. 

Matsumoto, M., Doi., M., 1994. Brownian dynamics simulation of DNA gel 
electrophoresis Mol. Sim. 12, 219-226. 

Noguchi, H., Takasu, M., 2001. Dynamics of DNA in entangled polymer solu- 
tions: An anisotropic friction model. J. Chem. Phys. 114, 7260-7266. 

Risken, H. 1989. The Fokker-Planck equation: Methods of solution and appli- 
cations. Springer Verlag, Berlin. 

Stellwagen, E., Stellwagen, N. C, 2002. The free solution mobility of DNA in 
Tris-acetate-EDTA buffers of different concentration, with and without 
added NaCl. Electrophoresis 23, 1935-1941. 

Stellwagen, E., Lu, Y., Stellwagen, N. C, 2003. Unified description of elec- 
trophoresis and diffusion for DNA and other polyions. Biochemistry 42, 
11745-11750. 

Streek, M. 2002. Migration of DNA on structured surfaces in an external field. 
Diploma thesis, Universitat Bielefeld. 

Tessier, F., Labrie, J., Slater, G. W., 2002. Electrophoretic separation of long 
polyelectrolytes in submolecular-size constrictions: A Monte Carlo study. 
Macromolecules 35, 4791-4800. 

Viovy, J.-L., 2000. Electrophoresis of DNA and other polyelectrolytes: Physical 
mechanisms. Rev. Mod. Phys. 72, 813-872. 



19 



